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Abstract 

In this paper, the modelling strategy of a Cosserat rod element (CRE) is addressed systematically for 
3-dimensional dynamical analysis of slender structures. We employ the exact nonlinear kinematic rela- 
tionships in the sense of Cosserat theory, and adopt the Bernoulli hypothesis. For the sake of simplicity, 
the Kirchoff constitutive relations are adopted to provide an adequate description of elastic properties 
in terms of a few elastic moduli. A deformed configuration of the rod is described by the displacement 
vector of the deformed centroid curves and an orthonormal moving frame, rigidly attached to the cross- 
section of the rod. The position of the moving frame relative to the inertial frame is specified by the 
rotation matrix, parametrized by a rotational vector. The approximate solutions of the nonlinear par- 
tial differential equations of motion in quasi-static sense are chosen as the shape functions with up to 
third order nonlinear terms of generic nodal displacements. Based on the Lagrangian constructed by the 
Cosserat kinetic energy and strain energy expressions, the principle of virtual work is employed to derive 
the ordinary differential equations of motion with third order nonlinear generic nodal displacements. A 
simple example is presented to illustrate the use of the formulation developed here to obtain the lower 
order nonlinear ordinary differential equations of motion of a given structure. The corresponding nonlin- 
ear dynamical responses of the structures have been presented through numerical simulations by Matlab 
software. 
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1 Introduction 



Three-dimensional slender structures undergoing large displacements and rotations are often encountered 
in various engineering systems such as vehicles, space structures, robotics, aircrafts, and microelectronic 
mechanical systems. Clearly, these systems consist of a set of interconnected components which may be 
rigid or deformable. For example, a typical MEMS device may consist of relatively heavy load bodies 
and thin springlike supports. For such a system, each heavy body can be assumed to be a rigid body 
and each springlike component can be described as a deformable body. Since each of interconnected 
components of such a system may undergo large displacements and/or rotations, an effective modelling 
strategy that addresses very well to the strongly nonlinear dynamic behavior is crucial in estimating 
system performance and guiding the reliability verification process. 

Nonlinear finite element method provides a general approach to structural modelling of multibody 
systems that consist of interconnected rigid and deformable components. A number of papers has re- 
cently been published, presenting new concepts and new algorithms for modelling highly flexible spatial 
frame structures [1, 2, 3]. An overview and comprehensive treatment of this topic can be found, for 
instance, in [4, 5]. The Cosserat approach, that can accommodate a good approximation the nonlinear 
behavior of complex structures composed of materials with different constitutive properties, variable ge- 
ometry and damping characteristics [6, 7, 8, 9], has been utilized to develop finite element formulations 
for deformable bodies. The finite element approach based on the Cosserat theory (geometrically exact 
finite-strain beam theory) is usually attributed to Reissner [10] and Simo [11]. Simo [11] has discussed 
a convenient parameterization of the rod model developed by Antman [12], and Simo and Vu-Quoc [13] 
have considered the associated finite element formulation. The computational procedure in [13] uses a 
variational formulation of the equations of motion and an expansion of the kinematic quantities in terms 
of shape functions and nodal values. Many modern finite element developers of the three-dimensional 
beam theories, e.g. Jelenic and Saje [14], Smolenski [15], and Zupan and Saje [16] based their approach 
on the geometrically exact beam theory. Another approach based on a system of Cosserat-type bodies can 
be traced back to the work of Wozniak [17]. Homogeneously deformable bodies have been analyzed as 
pseudo-rigid bodies [18] and Cosserat points [19]. The theory of a Cosserat point is a special continuum 
theory that models deformation of a small structure that is essentially a point surrounded by some small 
but finite region. The numerical procedure based on the theory of a Cosserat point proposed in [19, 20] 
has been used to study the dynamics of spherically symmetric problems in [21]. Recently, the theory of 
a Cosserat point has been generalized to model a fully nonlinear finite element for the numerical solution 
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of 3-D dynamic problems of elastic beams [22]. 

However, in practice the use of FEM codes to simulate complex multibody systems such as MEMS 
devices is prohibitively cumbersome, expensive, and time consuming. On the other hand, FEM models 
use numerous variables to describe the device state. This may lead the process of mapping the design 
space complex and the relationship between each of these variables and the overall device performance 
is not clear to designers. Recently, component level modelling methods, which contain a library of pa- 
rameterised behavioural models for frequently used MEMS components [23, 24], have been developed. 
In [23, 24], every component is described as a single element in contrast to FE models where the compo- 
nent is normally discretised into many elements. Consequently, lower degree models are established and 
the simulation time can be greatly reduced. The mechanical behaviour of the components, however, is 
often modelled using basic models, containing e.g. linear stiffness relationships and/or approximations 
of basic nonlinearities. 

Recently, motivated by the developments in MEMS modelling, the Cosserat theory has been employed 
to develop a novel modelling strategy that addresses very well to the practical needs for rapid modelling 
of slender structures such as the springlike components in MEMS, see Wang et al. [25]. This modelling 
strategy has been successfully used to investigate the non-ideal properties of typical MEMS beams [26]. 
In the sense of Cosserat theory, the motion of rods in three-dimensional space can be demonstrated by 
behaviors of a reference curve and three perpendicular unit vectors (directors). Consequently, the equa- 
tions of motion are nonlinear partial differential equations, which are functions of time and one space 
variable. For static problems, however, the equations become nonlinear ordinary differential equations, 
which can be solved approximately using standard techniques like the perturbation method to satisfy 
boundary conditions. In contrast, for dynamical problems, it is necessary to introduce a numerical pro- 
cedure which discretizes the equations. In the strategy for modelling of a Cosserat rod element [25], the 
basic kinematic quantities are the position of a point on the Cosserat curve and an orthogonal transfor- 
mation that define the rotation of an orthogonal triad attached to the cross-section at each point of the 
Cosserat curve. This enables description of a rod using nonlinear ordinary differential equations in terms 
of the generic nodal displacements of a CRE. 

As an initial consideration, the modelling strategy in [25] is developed for 2-D case. In this paper, 
the modelling strategy of CRE is addressed systematically for the 3-D problems. The fundamental prob- 
lem of any finite element formulation is the choice of the shape functions. The approximate solutions 
of the nonlinear equations of motion in quasi-static sense are chosen as the shape functions with up to 
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third order nonlinear terms of generic nodal displacements. In three dimensions, the nonlinear differen- 
tial equations cannot be integrated in a closed form even in the static sense, therefore the perturbation 
method is employed here to solve the system approximately. For the sake of simplicity, the Kirchoff con- 
stitutive relations are adopted to provide an adequate description of elastic properties in terms of a few 
elastic moduli. Based on the Lagrangian constructed by the Cosserat kinetic energy and strain energy 
expressions, the principle of virtual work is used to derive the ordinary differential equations of motion 
with third order nonlinear generic nodal displacements. The essential features and novel aspects of the 
present formulation for CREs are briefly summarized below. 

1. The shape functions for CREs are derived from the differential equations governing the flexural- 
flexural-torsional motion of extensional rods, taking into account all the geometric nonlinearities 
in the system. Consequently, the higher accuracy of the dynamic responses can be achieved by 
dividing the rod into a few elements which is much less than the traditional finite element methods 
in which the interpolation functions are usually extremely simple functions such as low order 
polynomials. 

2. The mathematical simplicity when formulating deformable bodies enables more convenient for 
modelling the multibody systems that consist of interconnected rigid and deformable components. 

3. The resulting nonlinear ordinary differential equations with lower degree-of-freedom are typically 
easy to simulate or integrate into system-level simulations. 

An outline of the main contents of this paper is as follows. We begin in section 2 by introducing 
the basic definitions and kinematic assumptions on the nonlinear elastic rods that can suffer flexure, 
extension, torsion, and shear. The rotational vector that is free both of singularities and constraints is 
employed as a parametrization to specify the deformed configuration space. We limited our attention 
here to the modelling of Cosserat rod elements in which the small effect of shear will be neglected. The 
governing equations of motion and the Kirchoff constitutive relations are presented in Section 3. The 
straightforward perturbation method is employed in Section 4 to solve the corresponding static problem. 
The approximate solutions obtained are subsequently used as shape functions of Cosserat rod elements. 
In section 5, Lagrangian approach is employed to formulate the nonlinear ordinary differential equations 
of motion of Cosserat rod elements. In terms of the shape functions derived in Section 4, the Lagrangian 
is constructed by the Cosserat kinetic energy and strain energy expressions, and the virtual work done 
by external point loads and distributed loads is discussed. A simple example is presented in section 
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7 to illustrate the use of the formulation developed here to obtain the lower order nonlinear ordinary 
differential equations of motion of a given structure. The corresponding nonlinear dynamical responses 
of the structure have been presented through numerical simulations by Matlab software. 

The following conventions and nomenclature will be used through out this paper. Vectors, which are 
elements of Euclidean 3-space M^, are denoted by lowercase, bold-face symbols, e.g., u, v; vector- 
valued functions are denoted by lowercase, italic, bold-face symbols, e.g., u, v; tensors are denoted 
by upper-case, bold-face symbols, e.g., I, J; matrices are denoted by upper-case, italic, bold-face sym- 
bols, e.g., M, K. The three vectors {ei,e2,e3} are assumed to form a fixed right-handed orthogonal 
basis. The summation convention for repeated indices is used. The symbols <9 t and d s denote differen- 
tiation with respect to time t and arc-length parameter s, respectively. The symbols (') and (' ) denote 
differentiation with respect to dimensionless time parameter r and dimensionless length parameter a, 
respectively. 

2 Kinematical Preliminaries 

2.1 Basic definitions and kinematic assumptions 

Adopt Cartesian coordinates (x, y, z) in inertial basis (ei, e2, 63) with Newtonian time t. According to 
the Bernoulli hypothesis the plane cross-sections suffer only rigid rotation during deformation and remain 
plane after deformation and preserve their shape and area. For the sake of convenience, we introduce the 
following definitions: (1) the reference configuration, where the geometrical and mechanical variables 
of the rod, including the loading, are known; (2) an arbitrary deformed configuration, where only the 
loading is prescribed, while the remaining variables are unknown. 

It is therefore convenient to introduce an orthonormal basis dj(s, t), (i = 1, 2, 3) of a cross-section 
at s, termed the moving basis, such that d3 is normal to the rotated cross-section, and di and d2 lie in 
the plane of the rotated cross-section. The motion of a rod segment can be modelled as a Cosserat rod 
whose configuration is described by its neutral axis r(s, t) (Cosserat curve) and 3 orthogonal unit vectors 
dj(s, t), (i = 1, 2, 3) (Cosserat directors) as shown in Figure 1. 

At any time, r describes the axis of the rod whose cross-section orientations are determined by d; such 
that d s r ■ d.3 > 0. This condition implies that (i) the local ration of deformed length to reference length 
of the axis cannot be reduced to zero since \d s r\ > 0, and (ii) a typical cross-section (s = sq) cannot 
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Figure 1 : A simple Cosserat model 



undergo a total shear in which the plane determined by di and d 2 is tangent to the curve r(-, t) at r(so, t) 
[10]. In general, as a result of shear deformations of the rod, the cross-sections are not perpendicular to 
the line of centroids. 

In an inertial Cartesian basis {ei, e2, 63} we may write 

r(s, t) = n(s, t) ei = x(s, t)e± + y(s, t)e 2 + z(s, t)e s . (1) 

The motion involves both the velocity of the curve, dtr(s, t), and angular velocity of the cross-sections 
w(s, t) so that 

d t di{s,t) = w(s,t) x di(s,t). (2) 

In a similar manner the strains of the Cosserat rod are classified into "linear strain" vector v(s, t) = 
d s r(s, t) and "angular strain" vector u(s, t) so that 

d s di(s, t) = u(s, t) x di(s, t). (3) 

It follows from the definition (2) that 

di x d t di = di x (w x d») = w(d, • di) - dj(dj ■ w) = 2w. 

Therefore, 

w = - di x d t di. (4) 
6 



Similarly, from the definition (3) we have 

u=-d;X<9 s d;. (5) 

Since the basis {di, d 2 , d%} is natural for the intrinsic description of deformation, we decompose rele- 
vant vector valued functions with respect to it: 

v(s, t) = Vi(s, t)di(s, t), u(s, t) = Ujdj(s, t), w(s, t) = Widi(s, t). (6) 



2.2 Parametrization of the rotation matrix 



There is a number of choices for the parametrization of rotation matrix, for example, the Euler angles, 
the quaternion parameters, and the rotational vector being the most usual [27]. Here, we employ the 
rotational vector that is free both of singularities and constraints. Because of the orthogonality the rota- 
tion matrix is a proper orthogonal matrix in SO(3), its nine components can be expressed by only three 
independent parameters. Denote S the spin matrix of a vector a = as 

-a 3 a 2 
5(a) = fl3 -ai • 

— 02 CL\ 

Then, the rotation matrix R is determined by the expression [27] 

k ?2/ 



(8) 



where = 0jej is the rotational vector, S(4>) is the spin matrix of defined by (7), and <fi = {<j>\ + 
<f>2 + ^l) 1 / 2 is the rotational norm or the length of the rotational vector. An expansion of trigonometric 
functions in Eq. (8) in MacLaurin's series yields 



R 



i + s + ±-s 2 + l-s 3 + --- + ±-s n + • 

2! 3! n! 



exp S. 



(9) 



Thus, the rotation matrix may alternatively be expressed by an exponential map, the exponentiation of 
the spin matrix associated with the rotational vector. Note that, as a consequence of the exponentiation of 
the spin matrix S(cf>) being equal to R(4>) G SO(3), the spin matrix S(<p) belongs to Lie algebra so(3) 
associated with the Lie group SO(3) [28]. 

Conversely, taking a given orthogonal matrix R as a rotation matrix, the associated rotation vector cj) 
can be derived from (7) and (8). The rotational norm cj> can be calculated by 

_ x Tr{R) - 1 



cos 



(10) 
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By taking the matrix logarithm of R we can obtain the skew-symmetric matrix S as following. 

S = logR=— f—(R-R T ). (11) 
2 sin <f> 

Therefore <j> = <Pi^i with cpi = -S23, 4>2 = Sys, and <f) 3 = -5i 2 . 

In terms of the rotational vector cf>, Eqs. (7) and (8) give the exact value of the current rotation matrix. 
Using truncated MacLaurin's series of various order in Eq. (9), approximate values of the rotation matrix 
are obtained and corresponding simplified theories can be derived. For example, a so called first order 
theory is obtained if small rotations are assumed so that the quadratic and higher order terms in Eq. (9) 
may be neglected. 



2.3 Specifications for the deformed configuration space 

The position vector r(s, t) denned by (1) is an element of Euclidean vector space g/ft . The orientation 
of the moving basis is represented by the rotation matrix, which is en element of the Lie group SO(3). 
Accordingly, the set of all possible configurations of the rod is defined by 

= {(r,fl)|r : s R : s -> SO(3)}. (12) 

This set is referred to as the deformed configuration space. The quantities r and R are termed the kine- 
matic quantities of the rod. Since the rotation matrix is related to the three parameters, the components 
of the rotational vector cj>, the Lie group SO(3) of rotation matrices is three-parametric, i.e. it may be 
viewed as being a 3-D nonlinear differentiable manifold. 

For a typical slender rod such as the components in MEMS, the effect of shearing deformation can 
be negligible, the cross-section of the rod is therefore assumed to be perpendicular to the tangent to the 
Cosserat curve, i.e. 

v(s,t) =d s r(s,t) = \d s r(s,t)\d 3 (s,t). (13) 

In this case, we write 

d 3 M) = |g ~ ^t^ 61 + ^(5,^2 + Ms,t)e 3 (14) 

with 

uf( S , t) + I/|( S| t) + Ufa, t) = 1. (15) 
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where v\ , v 2 and 1/3 can be written as 



d s x(s,t) d s y{s,t) 



and ^3(5, t) 



(16) 



|a s r( S ,t)|' |0 a r( a ,t)|' ~"~ |0,r(M)l 

by differentiating the position vector r(s, t) defined in (1) with respect to s. 

We assume that the directors {di, d2, (I3} can be obtained by the following way. First of all, we rotate 
directors {ei, 62,63} about e3 with an angle ip to obtain the directors {di, d2, 63}. Then, rotation 
matrix R a associated with the rotational vector cp a = ip e 3 can be written as 



R a = R(4> a ) 



cos ip — sin ip 
sin <p cos ip 
1 



(17) 



Next, we introduce a rotational vector 



4>b 



sin 1 \Jv( + 



v 2 di + 



sm 



v\ d 2 



which rotates the vectors {di, d2, 63} to {di, d2, d.3}. Here, we assume that vf + z/f 7^ 0. Other wise 
d3 = 63, this rotating procedure can be omitted. Let R b be the corresponding rotation matrix associated 
with the rotational vector cj>b- Then 



R b = R{cj> b ) 



ufu 3 +u% uiu 2 (u 3 -l) 



-Vl v 2 
Consequently, the moving directors are obtained as: 

(vfv 3 + z/f) cos <p v\v 2 {y 3 — l)sinyT 



Vl 

v 2 
»3 



(18) 



+ 



vf + vl 



ei 



+ 



[v\v z + v\ ) sin 99 t v\v 2 (v 3 - 1) cos ip 



vf + ^ 



+ 



Vf + 1/2 



e 2 - (v\ cos 92 + v 2 sm(p) e 3 , 



(19) 



do 



{vf v 3 + v%) sin ip i V\V 2 {y 3 — l)cos<£ 



z/ 2 + i/| 



+ 



i/ 2 + z/ 2 



ei 



+ 



(i/fi/3 + i/ 2 ) cos <p viv 2 (v 3 - 1) sin (p 



vf + vl 



vf + vl 



e 2 + (vi sin ip - v 2 cos tp) e 3 , 



vi&i + z/ 2 e 2 + 1/363- 



(20) 
(21) 
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Obviously (p(s, t) is a variable related to torsion of the rod. Expanding the directors in polynomials about 
vi,U2,(j) an d reserving the terms up to third order, we have 

di(s,t) « M. - ^<p 2 {s,t) - ^vf(s,t) - ^ui(s,t)u 2 (s,t)(p(s,t) \ e a 

+ ^¥>(s,t) - - ^z/|(s,t)^(s,t) - ^<p 3 (s,t)j e 2 

+ ^-^i(s,t) - 1/2(3, t)(p(s,t) + ^vi(s,t)<£ 2 (s,i)^ e 3 (22) 

d 2 (M) ~ f-<p(s,t) - ^i(s,t)u 2 (s,t) + ^i/f(s,t)^(s,t) + ^0 3 (s,t)^ ei 

+ f-i/a(«,t) +^(a,t)^(a,t) + ^z/ 2 (s,t)^ 2 (s,t)^ e 3 (23) 

d 3 (s,t)«z/ 1 (s,t)ei + ^ 2 (s,t)e 2 + ^l-^z/ 2 (s,t)-^ 2 2 (s,t)^ e 3 (24) 

For the sake of convenience to describe the displacements and rotations in the inertia frame, we regard 
directors dj(s,i) (i = 1,2,3) as those obtained by rotating inertial frame {ei,e 2 ,e 3 } with a rotation 
vector 

= (j) x (s, t)ei + ^(s, t)e 2 + 4> z (s, t)e 3 . (25) 

Now, based on the relations (19)— (21), utilizing the inverse procedure mentioned in Section 2.2, the 

rotational norm <p and the spin matrix associated with the rotation vector cj> in (25) can be derived from 

_! Tr{R b R a ) - 1 
(p = cos . (26) 

and 



S = \og{R b R a ) = -^—-(R h R a - R T a Rl). (27) 
2 sm <p 

Consequently, the approximate relations between ((j> x ,(f)y,(f) z ) and (vi } V2,<p), up to third order, are 
obtained as 

(/>x(s,t) = -u 2 (s,t) + ^cp(s,t)u 1 (s,t) - i ^ 2 (s,t) + vl(s,t) - ^ip 2 (s,t)J u 2 {s,t), 
<j) y (s, t) = u^s, t) + i ^(s, t)^ 2 (s, t) + i ^i/ 2 (s, t) + i/ 2 2 (s, t) - ^ ^ 2 (s, t)^ t), (28) 
^(s, t) = <p(s, i) - — (^ 2 (s, i) + i/| (s, i)) ip(s, t), 
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or equivalently, 

^(s,t) = 4> y {s,t) + ^(f> x (s,t)(f> z (s,t) - ^ {<Pl(s,t) + <p 2 y (s,t) + 4%(s,t)) <p y {s,t), 

< u 2 (s,t) = -0 x {s,t) + ^M s > t )Ms,t) + 1 (</> 2 x (s,t) + (f> 2 y (s,t) + cg{s,t)) Ms,t), (29) 

^ <p(s, t) = 4> z {s, t) + — {(f> 2 x (s, t) + 4> 2 y (s, t)) <t> z {s, t). 

These relations are very useful in solving the static problem and will be used below to derive the shape 
functions for CRD. 

3 The Governing Equations of Motion 

The dynamical evolution of the rod with density, p(s), and cross-section area, A(s) is governed by the 
Newton's dynamical laws: 

p{s)A(s)d tt r = d s n(s, t) + f (s, t), 
d t h(s, t) = d s m(s, t) + v(s, t) x n(s, t) + l(s, t), 

where 

n(s, t) = m(s, t)di(s, t), m(s, t) = m^s, t)dj(s, t) (31) 
are the contact force and contact torque densities, respectively; while 

h(s,t) = hi(s,t)di(s,t) (32) 

denotes the angular momentum densities; i(s,t) and l(s,t) denote the prescribed external force and 
torque densities, respectively. 

The simplest constitutive model is based on the Kirchoff constitutive relations which provide an ade- 
quate description of elastic properties in terms of a few elastic moduli. One may exploit the full versatil- 
ity of the Cosserat model by generating the Kirchoff constitutive relations to include viscoelasticity and 
other damping, curved reference states with memory and effects to prohibit total compression. 

The contact forces, contact torques and the angular momentum are given as 

n = K(v-d 3 ), m = J(u), h = I(w) (33) 
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where, according to the Kirchoff constitutive relations, the tensors K, J and I are described as 

K(s,t) = K ii (s,t){d i (s,t)®d i (s,t)), 

2 

J(s,t) = ^2 J ii (*,t)(d i (s,t)®d i (s,t)) + J 33 (s,t)(d 3 (s,t)®d 3 (s,t)), 

2 

I(s,t) = ^ I ij (s,t)(d i (s,t) ^ djis,^) + I 33 (s,t)(d 3 (s,t) d 3 (s,t)). 

The corresponding components are given as 

K 11 = K 22 = GA(s), K 33 = EA(s), 

J H = L(s) ^ dA > J 22 = L(s) E ¥ dA > 

J 33 = L{s) + ^ dA > J 12 = -J21 = L(s) E ^dA, 

7 H = L(s) dA ' J 22 = L(s) Pi 3 )? dA > 

> J 33 = J A (s) P( S )tt 2 + dA > J 12 = - J 21 = J A (s) P( S )CvdA, 

where E and G are the Young's modulus of elasticity and shear modulus respectively. 



(34) 



(35) 



4 Shape Functions for Cosserat Rod Elements 

For convenience, consider a uniform and initially straight rod element of constant length L, supported in 
an arbitrary manner at s = a = and s = b = L. It is assumed in the following that the static equilibrium 
of the rod corresponds to the situation where the directions of d3 and e3 are coincident with each other 
and di, d 2 are parallel to ei, e 2 , respectively. The principal axes are chosen to parallel ei, e 2 and e3. 
For the sake of simplicity, it will be assumed that the axes along the directors di, d 2 and d3 are chosen 
to be the principal axes of inertia of the cross section at s, and centered at the cross section's center of 
mass. Then, for a uniform rod with cross-section area A(s), we have J\ 2 = J 2 \ = 0, I\ 2 = I 2 \ = and 

' K n = K 22 = GA(s), K 33 = EA(s), 

Jn = E L( S ) V 2 dA, J 22 = E J A(s) i 2 dA, 

< hi=pJ A{ s)V 2 dA, I 22 = pf A(s) edA, (36) 
J 33 = G f A(s) (e + v 2 ) dA = f (J n + J 22 ), 

, J 33 = P L( S ) + V 2 )dA = I u + I 22 . 

Assume that the shape functions for a CRE satisfy the corresponding static equations of (30), i.e. 

d s ii{ S ) = (37) 
d s m(s) + v(s) x n(s) = 0, (38) 
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where the contact force and contact torque densities are 



(39) 



n(s) = m(s)di(s), m = K n vx, n 2 = K 2 2V 2 , 713 = ^33(^3-1), 
^ m(s) = mi(s)di(s), m 1 = J11U1, m 2 = J22U2, m 3 = J 33 n 3 . 

with u(s) = |dj(s) x i9 s dj(s), and dj(s) (z = 1, 2, 3) are given by (22)-(24). 

As mentioned in Section 2.3, for a typical slender rod as the components in MEMS, the effect of 
shearing deformation can be negligible, therefore the cross-section of rod is assumed to be perpendicular 
to the tangent to the Cosserat curve, i.e. the strain vector v(s, t) = |<9 s r(s)| d 3 (s) satisfies the form (13). 
Thus, v\ = v 2 = and tj 3 = |5 s r(s)|. Consequently, instead of m = K\\V\ and n 2 = K22V2, the 
contact forces n\ and ri2 follow from (38) 



—d s m 2 - 113ml + u\mz , d s mi - u s m 2 + U2m 3 
, and ri2 = 



i'3 



i'3 



(40) 



As a prelude to expanding the nonlinear shape functions to a form suitable for a perturbation analysis 
of the motion, it is useful to introduce some natural scales to obtain a dimensionless equation of motion. 
Introduce the dimensionless variables 



s r_ x V - z 

(T = —, r = —, x = —, y = —, z = —, t = u) Q t, 

Lq Lq Lq Lq Lq 



(41) 



where Lq and loq are the reference length and natural frequency yet to be determined later, respectively. 

Assume that the dimensionless generic nodal displacements (boundary displacements and rotations) 
at a = and a = L/Lq are 

i T 



eX a eY a eZ a e$ xa e<S> ya e<f> za 



(42) 



and 



Qb 



eX b eY b eZ b e$ xh e<S> yb e$ zb 



(43) 



respectively. Substituting (42) and (43) into (1), we obtain the boundary conditions for x, y and z as 



x(0) = eX a , y(0) = eY a , z(0) = eZ a , 
x(l) = eX b , y(l) = eY b , z(l) = I + eZ b , 



(44) 



where I = L/Lq is the dimensionless length of the rod element. Substituting (42) and (43) into (29), we 
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obtain the boundary conditions for v\, u 2 and <p as 




< 




(45) 



Treating e as a perturbation parameter which is the order of the amplitude of the displacement and can 
be used as a crutch in obtaining the approximate solution, the shape functions can be obtained by solving 
the static equations (37) and (38) with the corresponding boundary conditions (44) and (45) and also the 
restrictions (40) on the assumption of neglecting the effect of shearing deformation. To do this, we seek 
a straightforward expansion 



Substituting (46) into (37) and (38) associated with (40) and, because Xj, y~i, Zi and (pi are independent 
of e, set the coefficient of each power of e equal to zero. This leads to a set of linear ordinary differential 
equations which can be solved using the Frobeniu's method [29] under the corresponding boundary con- 
ditions (44) and (45). The solving procedure has been implemented in a MAPLE file [30]. Consequently, 
the approximate series solutions are obtained and the 1st order ones are 



< 



x{o) 
z(a) 



exi(a) + e 2 x 2 (a) + e 3 x 3 (a) H , 

eyi(o-) + e 2 y 2 (a) + e 3 m{a) + ■■■ , 
a + ezi(a) + e 2 z 2 (a) + e 3 z 3 (a) + • • 
e<pi(a) + e 2 ip 2 (a) + e 3 <£ 3 (cr) H . 



(46) 



£i(<t) = X a + $ ya a - (3X a - 3X b + 2l% a + l$ yb )-» 

+ (2X a -2X b + l% a + l% b )-p 
m {a) =Y a - <f> xa a - {3Y a - 3Y b - 2Z* SM - l$ xh ) ^ 

3 1 

+ (2Y a -2Y b -l$ xa -l$ xb )^- 
zi(<r) = Z a + (Z b - Z a )j 



(47) 



^i(a) = $ za + ($ 2fe - $ za ) T . 
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To investigate deflections up to 3rd order nonlinearity in e it is adequate to adopt the truncated (46) to e 3 
order terms. The high order terms (up to third order) which are polynomials of a, can be easily solved 
using a MAPLE programm [30]. For example, X2(cr) = Ci<r 5 + C2CT 4 + C3CT 3 + C^o 1 with 

Ci = ^£j^( Z b - Z a )(2X a - 2X b + l$ ya + l$ yb ). (48) 

Accordingly to the time-dependent, rod shape under the quasi-static condition is specified with the 
(slowly) time-varying nodal displacements and rotations. 



5 Equations of Motion for Cosserat Rod Elements 



In this section, the Lagrangian approach is employed to formulate the ordinary differential equations of 
motion of Cosserat rod elements. The generalized Hamilton's principle which, in its most general form, 
is given by the variational statement 



Jt! Jtl 



(49) 



where 3? is the total kinetic energy of the system, V is the potential energy of the system (including 
the strain energy and the potential energy of conservative external forces), S( ■ ) represents the virtual 
displacement (or variational) operator, and 5W is the virtual work done by nonconservative forces (in- 
cluding damping forces) and external forces not accounted for in Y. 

Assume that the time- varying dimensionless displacements at the ends (a = a/L Q and a = b/L ) of 
the element model are 

1 T 



and 



Qa{ T ) 



<lbl T ) 



X a (r) Y a (r) Z a (r) $ xa (r) $ ya (r) $ 2a (r) 



X b (r) Y b (r) Z b {r) <S> xb {r) $ yb (r) $ zfe (r) 



respectively. Then, the generalized displacement vector for the element can be described by 

i T 



Q e (r) 



(50) 



(51) 



(52) 



Consistent with the kinematic and constitutive assumptions described in Section 2 and Section 3 and 
the shape functions derived in Section 4, the kinetic energy per unit length is 



1 1 

7= -{ P Ad t r-d t r + I(w,w)} = - {pA^L 2 f • f + wgl(w, w)} 



(53) 
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where p and A are the density of rod and the area of cross-section of rod, respectively. According to (1) 
and (4), the velocity dtr(s, t), and the angular velocity of the cross-section can be derived as: 



and 



d t r = d t xei + d t ye 2 + d t ze 3 = uj L (xei + ye 2 + ze 3 ) = uj L r 



w = ^di x d t di = ^uj di x d; = w w, 



(54) 



(55) 



respectively. 

Under small strain conditions the strain energy per unit length of rod can be expressed in terms of the 
strain vectors u and v as: 

U = \ { J(U, U) + iY33(^3 " I) 2 } = \ {^2 J ( a > a ) + K ^3 " l)j (56) 



where the strain vector is 
1 



: -di x d s di = 7-j— dj x d- = -^-u and v 3 = \d s r\ = \r' 



V 3 . 



(57) 



Utilizing the time varying generic nodal displacements introduced in (50) and (51) instead of the 
static generic nodal displacements introduced in (42) and (43) respectively, the time varying generic 
displacements at any point within the element can be expressed as nonlinear functions of the length 
parameter a and the nodal displacement vector q e (r). Based on the nonlinear shape functions derived in 
Section 4, we have 

x = xi(cj, r) + x 2 (cr, r) + x 3 (a, r), 
V = m(<r,T) + 2/ 2 (ct, r) + £3(0-, r), 
z = a + zi(a,r) + z 2 ((t,t) + z 3 (a,r), 
_ <P = 0i{cr,T) + 02(o;T) + ip 3 (a,T). 
where the ith terms Xi,yi, Z\ and are ith order functions of the nodal displacement vector q e {r). For 
example, based on (47) the 1st order terms are 



(58) 



a 



x\ (cr, r) = X a (r) + $ yo (r)o- - (3X a (r) - 3X b (r) + 2l% a (r) + 1%^))^ 



(7 



(2X a (r) - 2X b (r) + l<f> ya (r) + l% b (r))-^ 



m(a,T) = Y a {r) - $ xa (T)a - (3F a (r) - 3F 6 (r) - 2Z$ xa (r) - /^(r))-^ 



(59) 



+ (2Y a (r) - 2Y b {r) - l$ xa ( T ) - l<S> xb (r)) 



a 



zi(<t,t) = Z B (r) + (Z b (r) - Z a (r)) y 
01 (a, r) = $, a (r) + ($ gb (r) - $*a(r))y 



(7 
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The high order terms (up to third order), as indicated in Section 4, can be easily obtained using a MAPLE 
program [30]. Consequently, the time varying generic displacements at any point within the element can 
be written as 



x = x(a,q e (T)), y = y(a,q e (T)), z = z(a,q e (r)), <p = <p(a, Q e (r)). 

This leads r = f (<r, q e (r)). Moreover, from (16), we have 

x'(a,q e {r)) y'(a, q e (r)) 

v\ = TT77 — rrrrr = i r ))> ^ = ft? — ff\t = ^ 0, <? t 

Substituting ^i(<7, q e (r)), 1*2 (f, q e ir)) and </?(cr, q e {r)) into the expressions (22)-(24) yields 

di = d i (<7,g e (r)), i = 1,2,3. 

Similarly, from (28), we have 

(px = 0x(cr,g e (r)), (f)y = (f)y(a,q e (T)), <p z = <f> g (a, q e (r)). 

It follows from (55), (57) and (62) that 

w = idj x di = w(ct, Q e (r)), u = ^d; x d- = u(a, q e (r)) 

Therefore, the kinetic energy density (53) and the potential energy density (56) are expressed as 

7 = 7(a,q e (r),q e (r)), U = U(a, q e (r)). 

Then, the Lagrangian defined in the classical form Jzf = — "V are obtained as 

if (q e , if) = 2T{cf, if) - Y(q e ) = f ( 7(a, q e , if) - U(a, q e ) ) L da. 

Jo 



(60) 



(61) 



(62) 



(63) 



(64) 



(65) 



(66) 



So far we have not precisely defined the type of loading. Let us assume that a load acting on the 
element is composed from three additive parts. The first one is the interaction of the neighbored elements. 
The second one is the external point (concentrated) loadings acting on the nodes. The last one represents 
a distributed load with fixed direction and prescribed intensity as mentioned in Section 2. In keeping 
with the load definitions in the principle of virtual work, the total load has to be defined with respect to 
the inertial basis because the generalized nodal displacements are defined with respect to them. Thus, let 
us denote 



(67) 





fL(r) 




' fUr) 




Itair) 




'Ur) 


fair) = 


flair) 


j fb — 


fyoir) 


J l a 


w 


1 l b — 






flair) 




Ur) 




%air) 




Ur) 
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be the interaction force and torque vector at nodes a = and a = I respectively. 
Similarly, the external point loadings are expressed as 



















f C a(r) = 




, ft(r) = 


fy b (r) 


r — 


w 


> l b — 


w 




Hair) 








l c za {r) 







while the distributed forces and torques may be expressed as 





r) 




4 


a, 


r) 




r) 


, v d = 




a, 


r) 










a, 





The virtual work done by the distributed load has the form 



d dxja, g e ) td dy{a, q e ) ed dz{a,q e ) 
^ Q qe -*-$y g qe 1-5, g qe 



d d<j) x (a, g e ) d d(py(a, q e ) d d<j) z (a, q e ) 



Lq da 



dq e 



dq e 



dq e 



5q e Lq da 



For the sake of convenience, let 



f e (r) 



fair) 

fl(r) 
1\{t) 



f c a (r) 
ft(r) 



f de (r,q e ) 



dx(a,q e ) d dMo, Q e ) \ 

+ --- + Vz\ T ) — I L oda. 



dq e 



Then, the total virtual work done by the three additive parts are 

5W = (f ie + f ce + f de f ■ 5q e . 



(68) 



(69) 



(70) 



(71) 



(72) 



(73) 



Substituting (66) and (73) into (49), taking variations using the chain rule, and integrating by parts, 
yield the generalized Lagrange equations of motion for the Cosserat rod element: 

d_ (dL \ _ dL_ 

dr \dqj J dqj 



/f(r) + /f(T) + /f(T,q e ) 



(74) 
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For a general configuration with nonzero generic nodal displacements q e , the ordinary differential 
equations of motion with up to third order nonlinearities of displacements and first order kinetic terms 
can be obtained as 

M e q e + K e q e + g e (q e ) = f e (r) + / ce (r) + f de (r, q e ), (75) 

where M e and K e are mass and (linear) stiffness matrices of the element model, g e (q e ) is a nonlinear 
vector with quadratic and cubic terms of q e . Since the mass of a typical rod, such as the springlike 
support component of MEMS, is very small comparing with the mass of the main device in practice only 
the first order kinetic terms are reserved in Equation (75). 

The detailed expressions of M e , K e and g e (q e ) have been implemented in a MAPLE program [30]. 
For the sake of illustration, the explicit expressions of M e , K e and g e (q e ) for a cantilever beam as a 
special Cosserat rod element are listed in Appendix. 

6 Dynamical Responses of rods by Cosserat Rod Elements 

6.1 Assembly of equations of motion for whole system 

We could analyze all of the types of systems consist of a set of interconnected components described in 
the introduction by using Cosserat rod elements for the deformable parts or subdivided members. Two- 
and three-dimensional frame structures require rotation-of-axes transformation for actions and displace- 
ments. For the sake of convenience, in this section we shall examine only the type of structure which 
is aligned with reference axes, using properties of the Cosserat rod element developed in the preceding 
sections. The analysis of the response of a number of complex structures is beyond the scope of this 
paper and will be presented in future publications. 

After stiffness, mass, and actual or equivalent nodal loads for individual Cosserat rod element are 
generated, we can assembly them to form the equations of motion for a whole system. We define global 
displacement vector q holding the displacement variables for all mesh nodes, such that 

q = [X x Yi Zi $,i $„i $ gl X 2 Y 2 Z 2 $ x2 % 2 $ 22 • • • ] T . (76) 

The equations of motion for the whole system can be constructed by simply adding the contributions 
from all the elements. In this way, expanding the matrix or operator for each individual element to make 
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them the same size as the system matrices or operators, we have 
and 



(77) 



i=i 



i=l 



i=l 



f C (r) = £ /V> «)=£ /' e ( r ' 9). (78) 

i=l i=l 

where n e is the number of elements. In Equation (77) M and K represent the system mass matrix 
and the system (linear) stiffness matrix. Similarly, the action vectors / c (r) and f d {r,q) are actual and 
equivalent nodal loads for the whole system. The contributions from the interaction forces and torques 
from all the elements must be of balance and the total action must be vanished. Then the undamped 
equations of motion for the assembled system become 



Mq + Kq + g(q) = / c (r) + f d (r, q). 



(79) 



This equation gives the system equations of motion all nodal displacements, regardless of whether they 
are free or restricted. 

In preparation for solving the nonlinear dynamic equations (79), as in the standard finite element 
procedure, we rearrange and partition it as follows 



Mff M fr 






M rf M„ 




q r 



+ 



Kff K fr 







+ 


9f(Qf,Qr) 




q T 







(80) 



f%T) + f d f (T,q f ,q r ) 
fr(^) + fr(r,q f ,q r ) 

in which the subscript / refers to free nodal displacements while the subscript r denotes restrained nodal 
displacements. If the support motions (at constraints) are zero, the equation (80) can be simplified to 



Mf f qf + Kffq f + g f (q f ) = f c f (r) + f d f{r,q f ), 



(81) 



and 



M rf q f + K rf q f + g r (q f ) = f c r (r) + f d r (r, q f ), (82) 
which can be used for solving the free displacements Qy(r) and support actions /°(r), respectively. 
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6.2 Simulation results and discussion for a simple cantilever 



A cantilever, as shown in Figure 2, is now presented as a simple example to demonstrate high accuracy 
and excellent performance of the proposed Cosserat rod elements. Numerical calculations based on (81) 
are carried out for a uniform horizontal cantilever of length L = 0.3m, of constant cross section with 
width B = 0.01m and thickness D = 0.005m. The mass density and the Young's modulus are assumed 
to be p = 3.0 x 10 3 kg/m 3 and E = 2.08 x 10 8 Pa. 
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p, E, A, L 




Figure 2: Schematic of a simple Cantilever 

Dividing the cantilever into n e elements of equal length, we can establish the nonlinear differential 
equations of motion (81) for solving the free displacements. In what follows, the natural frequencies 
of the linearized system are studied and used to compare with those derived from the classical beam 
theory presented in textbooks (see, for example [31]), and numerical simulations for the responses of the 
nonlinear dynamical system (81) under external harmonic excitations are performed with Matlab. 

Table 1 Flexural natural frequencies based on CRE approach and exact Continua method 





Flexural frequencies in ei-e3 plane 


Flexural frequencies in 


e2-e3 plane 


(rad/sec) 


CRE 


CBT 


| Error | (%) 


CRE 


CBT 


|Error| (%) 


1 


29.7607 


29.7665 


0.0197 


14.8827 


14.8833 


0.0036 


2 


186.358 


186.544 


0.0995 


93.2838 


93.2718 


0.0129 


3 


522.329 


522.329 


0.0000 


261.868 


261.164 


0.2692 


4 


1028.68 


1023.56 


0.5005 


516.914 


511.778 


1.0035 


5 


1707.74 


1692.01 


0.5155 


857.104 


846.007 


1.3118 



First, the flexural natural frequencies calculated in terms of the linearized equations of the nonlinear 
system (81) obtained by Cosserat element approach, together with the theoretical results obtained by 
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employing the classical beam theory (CBT) are given in Table 1. The flexural natural frequencies in 
both ei-e3 plane and e2-e3 plane, based on the CRE approach, showed their excellent convergency (the 
corresponding results listed in Table 1 are found, when only five Cosserat rod elements are used). 




Number of elements 

Figure 3: Convergency test for the first flexural frequencies in e2-e3 plane. -A-, u\\ -()-, t02', W3. 

Figure 3 represents the CRE convergency tests corresponding to the first three flexural natural fre- 
quencies in e2-e3 plane of the rod. As can be seen for the first frequency, the jerrorj is found to be very 
small (< 0.1%) even when only two elements are used. In fact the | error | for the first frequency is only 
0.4535% when just one element is used. For the second and third natural frequencies, the results are 
converging with approximately 0.1% error, when six elements are used. 

In the second part of this example, based on the derived nonlinear system (81), numerical simulations 
are performed to investigate the dynamic responses of the cantilever under harmonic excitations. The dif- 
ferential equations of motion are full coupled by the nonlinear terms and could exhibit internal resonance 
introduced by the nonlinearities. They also exhibit external resonances when the external excitation is 
periodic and the frequency of a component of its Fourier series is near one of the natural frequencies of 
the system, or near a multiple of a natural frequencies. The detailed analysis of complex dynamic behav- 
ior, such as bifurcation and chaos, of the system is not the main focus of this paper. We only compare 
here the responses of the system, when different number of elements are used. 

The displacement and angular time histories of the free end of the cantilever under external loads 
fx(t) = 0.01cos(8t), fy(t) = 0.005 sin(8t) and at zero initial conditions are shown in Figure 4 when 
two elements are used and in Figure 5 when ten elements are used, respectively. It is interesting to 
note that amplitudes and periods of the responses are very closed in this two situations. To enhance 
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4 6 
t [sec] 



Figure 4: Displacement time histories of the rod with external loads f x (t) = 0.01 cos (8 * t), f y (t) 
0.005 sin(8 * t) and zero initial conditions: two elements case. 




4 6 
t [sec] 



Figure 5: Displacement time histories of the rod with external loads f x (t) = 0.01 cos(8 * t), f y (t) = 
0.005 sin(8 * t) and zero initial conditions: ten elements case. 

this observation, the phase plane diagrams for Y(t)-Y(t) in four different cases, namely one element, 
two elements, three elements and ten elements, are plotted in Figure 6 (a), (b), (c) and (d), respectively. 
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Comparing the four diagrams in Figure 6 shows that the modal when two or three elements are used can 
exhibit almost the same behavior of the model when ten elements are used. 
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Figure 6: Phase plane diagram of Y-Y with external loads f x (t) = 0.01 cos(8*£), f y (t) = 0.005 sin ( 
t) and zero initial conditions. 



According to the analysis of natural frequencies and the analysis of harmonic responses of the estab- 
lished nonlinear dynamic systems, we believe, in practical engineering problem, especial for the structure 
composed of springlike flexural components such as the device in MEMS, only a few Cosserat rod ele- 
ments are needed to model a flexural component. For the very slender flexural components, we can even 
use only one element to model such a component. 



7 Conclusion 

A Cosserat rod element formulation for the modelling of 3-dimensional dynamics of slender structures 
has been proposed in this paper. The modelling strategy of this new approach employed the exact nonlin- 
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ear kinematic relationships in the sense of Cosserat theory, and adopted the Bernoulli hypothesis. Finite 
displacements and rotations as well as finite extensional, torsional, and bending strains are accounted 
for. The Kirchoff constitutive relations, which provide an adequate description of elastic properties in 
terms of a few elastic moduli, are adopted. A deformed configuration of the rod is described by the 
displacement vector of the deformed centroid curves and an orthonormal moving frame, rigidly attached 
to the cross-section of the rod. The position of the moving frame relative to the inertial frame is specified 
by the rotation matrix, parametrized by a rotational vector. The approximation solutions of the nonlinear 
partial differential equations of motion in quasi-static sense are chosen as the shape functions with up to 
third order nonlinear terms of generic nodal displacements. This lends the approach very well to achieve 
higher accuracy of the dynamic responses of the model by dividing the slender rod into a few elements. 
Based on the Lagrangian constructed by the Cosserat kinetic energy and strain energy expressions, the 
principle of virtual work is employed to derive the ordinary differential equations of motion with third 
order nonlinear generic nodal displacements. 

A cantilever as a simple example has been presented to illustrate the use of the formulation developed 
here to obtain the lower order nonlinear ordinary differential equations of motion of a given structure. 
The natural frequency analysis for the linearized equations and the numerical simulation analysis for 
the nonlinear model show that in practical engineering problem, especial for the structure composed 
of springlike flexural components such as the device in MEMS, only a few Cosserat rod elements are 
needed to model a flexural component. 

The mathematical simplicity when formulating deformable components enables more convenient for 
modelling the multibody systems that consist of interconnected rigid and deformable components. The 
Cosserat rod element approach therefore is feasible to be used to capture the most significant character- 
istics of a multi- rigid and deformable body system in a few variables governed by nonlinear ordinary 
differential equations of motion. 

As the first step to present the Cosserat rod element approach, we have limited our attention to the 
modelling of Cosserat rod elements in which the effect of shear has been neglected. The extension 
of the present formulation to the modelling of more general Cosserat rod elements in which the finite 
extensional, torsional, bending strains as well as shear are accounted for is highly desirable. 

Acknowledgements The authors are grateful to the EPSRC (Computational Engineering Mathematics 
Programme) and the EC(Framework Programme) for financial support in this study. 
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Appendix 



Let us assume that a uniform cantilever beam of length L, of constant cross section with area A and 
density p, is fixed at s = and free at s = L. In this case, we have q a = 0, thus q e = q b . Consequently, 
M e , K e become 6x6 matrices, and g e (q e ) is a six dimensional nonlinear vectorial functions of q e = q b . 
They are 
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and 



91 {(f) = gi,\X b Z b + g 1>2 Y b ^ zb + 51,3^6*2/6 + 01,4**6**6 + Sl,5*6 + QlfiXfeyt 

+ gi, 7 X b Y b 2 + cft^Y^ + 5i, 9 X 6 Z 6 2 + si.w-X'ftfcirf, + ffi,iiX 6 ^ 6 

+ 51,12X5$^ + g hn Y b 2 <5> yb + 51,14^6^$^ + gi,i 5 Y b <f> xb <f> yb + g hW Z^ yb 

+ gi A7 Z b <$> xb <$> zb + gi tl8 $l b <S> zb + gi^ yb + gi t20 $ yb <!> 2 zb , 



(85) 



g2{q e ) = g2,iX b ® zb + g 2 , 2 Y b z b + g2,*z b § xb + g 2 ^ yb ^ zh + g2,bX b 2 Y b + g 2 , G x^ xb , 
+ g 2 jX b Y b <Z> yb + g 2j8 X b Z b <f> zb + g 2 , 9 X b <f> xb <f> yb + g 2 ,i Y b 3 + g 2 ,xiY b 2 <S> xb 

+ ff2,12>6^6 2 + 32,13^6^6 + #2,14*6^6 + 52,15^6*^6 + 52,16^fe*xfe 



+ g 2 ,i 7 Z b <5> yb <5> zb + g 2 ,i 8 % b + g 2 ,i9$ xb % b + 0j,2o*x&**& , 



(86) 
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53 (g 6 ) = 53,1^6 + 53,2^6*2,6 + 53,3^ + 53,4*6**6 + 53,5**6 + 53,6^6 + 93,7*% Z b 

+ g 3fi X b Y b <& zb + g 3t gX b Z b $ yb + gz,wX h $ xb $ zb + g 3A1 Y b Z b + g 3t i 2 Y b Z b $ xb 

+ 53,13*6*2/6**6 + 93,uZ b ^ 2 xb + g 3 ,15Z b $l b + 53,16**6*2/6**6 , (87) 

54 (<? 6 ) = 54,1^6**6 + 54,2*6-^6 + 54,3-^6**6 + 54,4*2/6**6 + 54,5^6 *6 + 54,6^6 *a;6 

+ g4jX b Y b $ yb + g4fiX b Z b $ zb + 54,9^6**6*2,6 + 04,io*b 3 + 54,ll5"f*r6 

+ 04,12*6^6 + 54,13*6**6 + 54,14*&* 2 ,6 + 54,15*&* 2 6 + 54,16^6 *x6 

+ 54,17-^6*2/6**6 + 54,18**6 + 54,19**&* 2 y& + 54,20**6**6 . ( 88 ) 

55 (<f) = 51,1^6^6 + 55,2*6**6 + 55,3^6*2/6 + 55,4**6* 2 6 + 55,5^6 + 55,6^6 *2/6 

+ g 5J X b Y b 2 + 95,8^x6 + 55,9^6^ + ffl,l ^6*x6 + 51,11^6*^6 

+ 55,12^6* 2 6 + 55, 13y 6 2 $ yb + 55,14*6-^6**6 + 55,15*6**6*2/6 + 55,16^*2/6 

+ 55,17^6**6**6 + 55,18**6**6 + 95,19*?6 + 55,20*2/6*^ , (89) 

56 (<? £ ) = 56,1^6*6 + 56,2^6**6 + 56,3*6*2,6 + 56,4**6*?/6 + 56,5^ 2 **6 + g6,eX b Y b Z b 

+ g 6J X b Z b ® xb + 56,8^6*2/6*^6 + 56,9*f**6 + 56,10^6^6*2/6 + 56,11*6*Z&*^ 

+ 56,12^6**6*2,6 + 56,13**6**6 + 56,14* 2 ,6* z 6 , (90) 

where, the coefficients of second order nonlinear terms are specified as: 

6(K 33 l 2 - 20J 22 ) 6(J 22 - Jn) 

51.1 — 2 53,i — ^4 , 51,2 — 52,1 — 56,1 — p 1 

K 33 l 2 - 60J 22 4J n - J 22 - J 33 

51.3 — 53,2 — -55,1 — ^3 j 51,4 — 54,1 — "56,2 — p , 

6(K 33 Z 2 -20J n ) K 33 l 2 -WJn 

52.2 = ^53,3 = ^4 , 52,3 = 53,4 = 54,2 = ^3 , 

Jn — 4J 22 + J33 1 1 _ -^33 

52.4 — —55,2 — 56,3 — J2 j 53,5 — 53,6 — 2 54,3 — ~ 2 55,3 



54,4 — —55,4 — 56,4 



j2 ' 30,0 30,0 2^^'° 2^°!° 

Jn — J22 



I 

The coefficients of third order nonlinear terms are specified as : 

18(7i^| 3 / 4 - 160J 22 K 33 l 2 - 560Jf 2 ) 
9l ' 5 ~ 175K 33 F ' 

_ 9(7K 2 3 l A - 260J 22 K 33 l 2 - 3360Jf 2 ) 
9lfi ~ 350K 33 / 6 ' 

= 18(7if 33 / 2 - 80(J n + J 22 )) 18(10K 33 / 2 (Jn - J 22 ) 2 + 112J n J 22 J 33 ) 
91 > 7 175/5 35J 33 K 33 F 

3(7K 33 Z 2 - 480J n + 220J 22 ) 18(10K 33 / 2 ( J n - J 22 ) 2 + 112J n J 22 J 33 ) 
9l ' 8 1T5Z 4 35J 33 K 33 /6 
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iff 3 / 4 + 840J 22 if 33 - 25200J| 2 
700 J 22 / 5 ' 

UK 33 l 2 - 500 J u - 80J 22 + 175 J33 52K 33 l 2 (J u - J 22 ) 2 + 504J n J 22 J 33 
175P 35J 33 K 33 l^ 

63JCf 3 / 4 - 520J 22 JC 33 / 2 - 38640J| 2 
700^ 33 Z 5 ' 

20 J n - 16Jn J 22 ~ 4Jn J33 ~ 4J| 2 + 4J 22 J 33 ~ J33 

5Jii/ 3 

3(7i^ 2 - 480J 22 + 220J n ) 9(10ir 3 3*Vii - J22) 2 + 112J n J 22 J33) 
350/ 4 35J33K33/ 6 

12( Jn - J 22 ) 
/ 4 

7if 33 Z 2 + 900(Jii + J 22 ) - 700J 33 118/f 33 / 2 (J n - J 22 ) 2 + 1428Jn J22J33 
700/ 3 + S5J 33 K 33 l 5 

Kj 3 l 4 - 8400 J| 2 
1400 J 22 / 4 ' 

BJnKssf - 2J 22 K 33 l 2 + J 33 K 33 l 2 - 240 J n + 60 J n J 22 + 60 Jn J33 

60Jn/3 

7if 33 / 2 - 240J n - 30J 22 40if 33 / 2 (J n - J 22 ) 2 + 462J n J 22 J 33 
1050/ 2 35J 33 Js: 33 / 4 

7if 33 / 4 - 270J 22 K 33 l 2 - 13860Jf 2 
1050K 33 / 4 

10 J n - 16Jn J 22 + Jn J33 - 4J| 2 + 4J 22 J 33 - J| 3 ^ 
10J n / 2 

6(7Jf 33 / 2 - 480J 22 + 220Jn) 18(10if 33 / 2 ( Jn - J22) 2 + 112 Jn J22 J33) 
350/ 4 35J 33 K 33 l e 

18(7iff 3 / 4 - 160JiiJf 33 / 2 - 560 J n ) 

175K 33 / 7 ' 
9(7Kf 3 Z 4 - 26OJ11K33/ 2 - 3360J 2 !) 



35OK33/ 6 
Ky A + 840JiiK 33 - 25200 J n 

700 J u l 5 ' 
63A'f 3 / 4 - 520JnA 33 / 2 - 38640J n 

700A 33 1 5 ' 
14A 33 / 2 - 500 J 22 - 80 J n + 175 J 33 52K 33 l 2 (J u - J 22 ) 2 + 504J n J 22 J 33 



175J 33 JC 33 / 5 35J 33 A 33 / 5 
22 — 16Jn J 22 — 4J 22 J 33 — 4J n + 4Jn J 33 — J 33 



5J 22 / 3 

Kl 3 l 4 - 8400 J n 
1400 Jn/ 4 
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52,17 
52,18 
52,19 
52,20 
53,14 

53,15 
53,16 
54,18 
54,19 
54,18 
55,19 
55,20 

and 



hJ 22 K 33 l 2 - 2J u K 33 l 2 + J 33 K 33 l 2 - 240Jf 2 + 60Jn J 22 + 60J 22 J 33 



60J 22 / 3 



7K 33 / 4 - 270JnK 33 / 2 - 138 60 J n 



lOSO^/ 4 

7if 33 / 2 ~ 240J 22 - 30Jn 40K 3 3^ 2 (Jii - J22) 2 + 462Jn J 22 J 33 

1050Z 2 35Z 2 
10J| 2 - 16Jn J 22 + J 22 J 33 - + 4Jn J 33 - J| 3 
10J 22 / 2 

if 33 (ll^ 2 ~ 840 J n ) 
6300 Jul 

K 33 (llK 33 l 2 - 840 J 22 ) 
6300 J 22 / 

-^33(2^11 — J11 J33 ~ 2Jf 2 + J22J33) 

120JnJ 22 
7if 33 / 4 - I8OJHK33/ 2 - 7560 J n 

1575^33^ 3 ' 
UK 33 l 2 - 180(J n + J 22 ) + 175J 33 2 85K 33 l 2 (J n - J 22 ) 2 + 30 24J n J 22 J33 

1575/ 315J 33 J^33/ 3 
12J 2 X + 28Jn J 22 — 12Jn J33 — 20Jf 2 + 2J 22 J33 + 3J| 3 

60J 22 / 

7i^ 4 - 180J 22 K 3 3^ - 7560J 2 2 



1575^33^ 

12J| 2 + 28Jn J 22 — 12J 22 J33 — 207^ + 2Jn J33 + 3J| 3 



60J U Z 



52,5 = 


= 51,7, 


52,6 = 


551,8, 


52,8 = 


= 51,14, 


52,9 = 


= 51,15, 


53,7 = 


= 51,9, 


53,8 = 


= 51,14, 


53,9 = 


551,16, 


53,10 


= 51,17, 


53,11 


= 52,12, 


53,12 


= 251,16, 


53,13 


= -52,17, 


54,5 = 


551,8, 


54,6 = 


= 5i,io, 


54,7 = 


= 51,15, 


54,8 = 


= 51,17, 


54,9 = 


= 351,18, 


54,10 = 


r —|52,ll, 


54,11 


= 52,13, 


54,12 


= 52,16, 


54,13 


= 352,18, 


54,14 


= 52,19, 


54,15 = 


= 52,20, 


54,16 


= 53,14, 


54,17 


= 53,16, 


55,5 = 


= §01,6, 


55,6 = 


= -51,11, 


55,7 = 


-552,7, 


55,8 = 


= —51,15, 


55,9 " 


= —51,16, 


55,10 


= -51,18 


55,11 


= -35i,i9, 


55,12 = 


r -51,20, 


55,13 


= -52,14, 


55,14 


= —52,17, 


55,15 


= —52,19 


55,16 


= 53,15, 


55,17 = 


r -53,16, 


55,18 


= —54,19, 


56,5 = 


= 51,12, 


56,6 = 


= 51,14, 


56,7 = 


= 51,17, 


56,8 = 


251,20, 


56,9 = 


= 52,15, 


56,10 


= —53,13, 


56,11 


= 252,20, 


56,12 


= 53,16, 


56,13 = 


= —54,20, 


56,14 


= 55,20- 
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